Attractive forces between anisotropic inclusions in the membrane of a vesicle 
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The fluctuation-induced interaction between two rod-like, rigid inclusions in a fluid vesicle is 
studied by means of canonical ensemble Monte Carlo simulations. The vesicle membrane is repre- 
sented by a triangulated network of hard spheres. Five rigidly connected hard spheres form rod-like 
inclusions that can leap between sites of the triangular network. Their effective interaction poten- 
tial is computed as a function of mutual distance and angle of the inclusions. On account of the 
hard-core potential among these, the nature of the potential is purely entropic. Special precaution is 
taken to reduce lattice artifacts and the influence of flnite-size effects due to the spherical geometry. 
Our results show that the effective potential is attractive and short-range compared with the rod 
length L. Its well depth is of the order of k/10, where k is the bending modulus. 

PACS: 87.22.Bt,68.35.Md,02.70.Lq 



I. INTRODUCTION 



Lipid membranes are interesting systems in statistical 
physics and are the subject of many theoretical and ex- 
perimental investigations because of the abundance of ef- 
fects they exhibit. Often fluid membranes are considered, 
which feature the unusual combination of finite bend- 
ing stiffness and vanishing in-plane shear stress. This is 
caused by the constituent lipid molecules which can be 
sheared against each other, but resist bending normal to 
the plane Q] . 

Biological membranes also contain inclusions, i. e. im- 
purities that differ from lipid molecules chemically and 
mechanically. Inclusions are embedded in the mem- 
brane and can diffuse laterally. Examples range from 
rather large inclusions such as proteins and polymers to 
very small bodies such as so-called gemini comprising 
two lipid molecules whose head groups are chemically 
bonded M. In general, any membrane component that 
deviates in its mechanical properties from lipid molecules 
will be called inclusion. Most of these are rather rigid 
and thus, the presence of an inclusion locally stiffens the 
ambient membrane. Inclusions diffuse within the mem- 
brane with typical speeds of a few microns per second M . 

Forces between membrane inclusions currently receive 
considerable interest P-pj. They can be divided into 
two classes, direct forces due to electrostatic and van der 
Waals interactions, and indirect forces which are medi- 
ated by membrane fluctuations. The latter are of inter- 



est here. However, indirect fluctuation forces between 
membrane inclusions should not be confused with deple- 
tion forces [mO| (although both are entropic in origin) 
that exist when small particles are depleted from the gap 
between bigger ones. Indirect inclusion interactions were 
investigated theoretically, both for isotropic (rotationally 
invariant in the membrane plane) p-p 11,13], anisotropic 
(e. g. rod-like) inclusions MM, ^^nd under lateral mem- 
brane tension |^ . The above quoted contributions focus 
on the range ^ 3> s 3> i, where s is the center-center 
distance between two inclusions whose linear, in-plane 
size is L, see Figure || and ^ is the persistence length, 
that is the distance in the membrane over which the cor- 
relation of the surface normals decays |Q. Based upon 
perturbative approaches, it was found that there is an 
(attractive) long-range interaction potential of the form 
$ ex —ksT (L/s)'^, both for isotropic M and anisotropic 
y,|| inclusions (fee and T are Boltzmann's constant and 
temperature, respectively). However, its magnitude was 
predicted to be much smaller than /cbT over the range of 
s:^ L. 

Short range attractive interactions (s < L) have been 
predicted by Netz who treated the interactions be- 
tween stiff inclusions in a membrane analytically. For his 
model Netz obtains a logarithmically decaying, attractive 
interaction potential whose magnitude increases with in- 
creasing membrane stiffness (see Sec. ||,0). Other stud- 
ies dealing with short-range interactions between mem- 
brane inclusions were performed by Dan et al. |13|| and 
Aranda-Espinoza et al. ||l3] who considered membrane 
Hamiltonians with contributions from compression (ex- 
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pansion), spontaneous curvature of the membrane, and 
bending stiffness. These calculations have been carried 
out in the limit of vanishing temperature where mem- 
brane fluctuations do no longer exist. Also, s in ||l^,|ljl is 
typically of the order of the membrane thickness, whereas 
the papers ^-^ assume a thickness much smaller than 
s. 

The present paper is also concerned with short-range 
interactions but between rod-like inclusions for s ~ L 
embedded in the surface of a vesicle. Throughout this 
paper we are exclusively concerned with the case of zero 
spontaneous curvature. To the best of our knowledge, 
numerical simulations for three-dimensional fluctuating 
membranes with finite bending stiffness and inclusions, 
which are considered here, have not yet been carried out. 

The remainder of this paper is organized as follows. 
Model and simulation algorithm are detailed in Sec. ||. 
In Sec. P|, pair distribution functions are introduced. 
Details of their numerical determination are presented in 
Sec. [V. Sec. M is devoted to a presentation of the results 
obtained in this work. The paper concludes in Sec. VI 
with a summary and discussion. 




FIG. 1. Curved membrane with two rod-like surface in- 
clusions of length L and width e <^ L, separated by the cen- 
ter-to-center distance s. The rods are rotated by the in-plane 
angles 9i , 62 against the connection vector. Sketches a) and b) 
show two rods in parallel settings. A perpendicular arrange- 
ment similar to c) will be called "T- formation" , a setting such 
as d) "in-line formation". 



II. THE MODEL 



The approximations employed in ^-^ are only valid 
in the range s ^ L and break down for the interesting 
case oi s ~ L. Also, the membrane can no longer be 
regarded as a continuous surface, as the discrete lipid 
network becomes more and more influential. Hence, nu- 
merical simulations are required to obtain quantitative 
results. In this regard the Monte Carlo method provides 



a particularly powerful technique by which thermophys- 
ical properties of equilibrium systems can be computed 
in a rather simple and straightforward manner |1^ . 

The model membrane consists of N hard spheres of 
diameter a, connected by rigid bonds (tethers) of length 
r with a < r < \/3a. This model has been employed pre- 
viously |lq-pj|]. Membrane fluidity is realized by the 
bond-flip algorithm first proposed in ||20| . On a closed 
triangular network, each bond can be regarded as one 
of the diagonals in the quadrilateral formed by the four 
surrounding bonds. The bond-fiip algorithm rotates the 
bond within this quadrilateral, so that it represents the 
other diagonal after the operation. This method allows 
for vertex diffusion, as any triangulation can be trans- 
formed into any other pC| | . The bending energy is com- 
puted from the Helfrich Hamiltonian pl| by integrating 
over the surface S 



n= lds{--H^ 



RK 



(1) 



where i/ = ci + C2 is the sum of the principal curva- 
tures on the surface, K ^ ci C2 the Gaussian curvature, 
and K, K, are, respectively, the bending modulus and the 
Gaussian modulus. The term /d5 K is constant in fixed 
surface topology due to the Gaufi-Bonnet theorem pi[ . 
The term /d5 H^ is discretized as in [Q , 



dS* iJ^ « 47r K 



Nb 



i=l 



1-nA HA, 



(2) 



where the sum runs over all Nb = 3(iV — 2) bonds in the 
network. Each bond has two adjacent triangles Aa{i) and 
At,(i), whose outer unit normals are denoted by ha^ and 




FIG. 2. Portion of a triangular network of hard spheres 
with two inclusions, consisting of five vertices each. To 
the ambient membrane, vertices and bonds within inclusions 
appear "frozen" and cannot be moved (or flipped) by an 
MC-step. An inclusion can leap as an entity (arrow) by re- 
leasing the old vertices and freezing a new set (bold circles), 
after having moved the vertices. Pi is the position of the 
selected vertex from which the new rod is tried to be placed 
along the tangential direction u, see text. 



The simulation algorithm consists of two independent 
parts. The first one follows the one described in ^^ . One 
Monte Carlo step (MCS) consists of attempting to move 
N randomly selected vertices to new positions within a 
cube [—p,p]'^, centered at their current positions. Next, 
N bonds are selected and attempted to be flipped. Both 
processes are accepted or rejected on the basis of their 
associated change in energy according to the Metropolis 
algorithm |23]. 

In addition, mobile, rod-like inclusions are embedded 
in the membrane. They consist of five rigidly connected 
hard spheres. None of the spheres in this "frozen" rod- 
like configuration can be moved independently, nor can 
the four inner bonds be fiipped. A rod and its ambi- 
ent membrane are thus similar to a spine in which some 
of the vertebrates are fused. Their biological equivalent 
would come closest to (multi-)gemini, since the inclusion 
constituents are equal to those of the plain membrane. 

The inclusions can move as an entity, thereby resem- 
bling lateral diffusion. This is accomplished by the "rod- 
leap" algorithm, described below, which is carried out 
in the second step of the simulation. A new set of five 
vertices is selected, which may overlap with the old one. 
To find the new vertices, a random vertex next to the old 
rod is selected and the four remaining vertices are chosen 
along the direction of a random unit tangential vector u. 
Subsequent vertices in the rod must be connected by a 
bond, as shown in Figure ||. Then a Monte Carlo move 
is attempted to simultaneously move the last four of the 
new vertices to the new positions 



Pi + (u-P„)u, 



(3) 



where Pi = Pi is the position of the selected vertex and 
the P„, {n = 1 . . . 5) are the old vertex positions. The 
orientation of the local tangential plane next to Pi is 
found by averaging all the neighboring triangle normals. 

The effect of placing a rod somewhere is to flatten 
the membrane locally. Note that this does not affect 
the surface topology. No holes or contact angles at the 
membrane-inclusion boundary are introduced and hence 
the Gaufi-Bonnet theorem remains valid. If finite con- 
tact angles are assumed |§,|2| , the number of the inclu- 
sions must be kept constant, whereas in the present case 
the GauB-Bonnet theorem would remain valid even if the 
number of inclusions were varied. This flattening effect 
quenches fluctuations normal to surface of the vesicle. 
On the other hand, inclusions hardly increase the net 
bending energy of the ground state vesicle (sphere) of 
Eo = Snu [|l^. Figure g shows a snapshot of a simula- 
tion with N = 1012 vertices. 

The total number of vesicle vertices iV and number of 
rods Nft must be carefully chosen to ensure that the vesi- 
cle is not too strongly perturbed by the rods. Also, the 
impact of multi-body effects among the rods is crucial for 
the reliability of the simulations. The rods on the mem- 
brane can be considered a lattice gas that should be as 
dilute as possible. In the present simulation, Nji = 4 



rods were placed on a vesicle of TV = 1012 vertices. 
With this choice, only 2% of the vertices are occupied 
by inclusions. The average rod length of i — 4(r), 
where (r) « a(l + \/3)/2 is the mean bond length, is 
only about half as long as the average vesicle radius 
of gyration that is computed to be {Rg) = 8.13 (r) for 
n/k^T = 7.5. An alternative measure of coverage of the 
membrane by inclusions is obtained by assigning a disk 
of size Al = TT (i/2)^ to each rod. The total area cov- 
ered by the disks NrAl/ {A) is about 6% of the average 
vesicle surface area {{A) k, 841 (r) ). Thus, we conclude 
that with the above choice of constants, the vesicle is 
not strongly perturbed by the inclusions and many-body 
effects are expected to be negligible (see also below). 




FIG. 3. Simulation snapshot of a vesicle with A^ = 1012 
vertices and K/ksT = 7.5. The shades on the surface may 
be thought of as reflections from a spotlight. Two inclusions 
with their incorporated spheres are visible on the membrane. 

Thermal fluctuations induce a finite effective surface 
tension r in fluid membranes. Therefore, in principle r 
introduces an additional length scale. However, this ten- 
sion was estimated as r « k^T^^"^ ||2^ and thus is very 
small in the regime that is of interest here. In the case 
of biological membranes, experiments show [gj| that the 
effective surface tension is negligibly small. 



III. PAIR DISTRIBUTION FUNCTIONS 

The goal of the present work is the computation of the 
fluctuation-induced interaction potential of the inclusions 
'1'rr(sj ^ii ^2, i"^)- The total energy of a given triangula- 
tion TJv, as given by ([I]) and (|), depends on the positions 
of the vertices Xi . . . Xjv and on the N x N connectivity 



(adjacency) matrix S, where S^ = 1 if vertices i and j are 
connected by a bond and zero otherwise. The partition 
function can then be written as |ial 



dTN exp{-f3{n[TN] + C/bond)) 



N 



IV. NUMERICAL DETAILS 

We approximate g^^^ by a histogram g^- . To do so, a 
rod pair histogram riij of 25 x 25 entries is generated. Af- 
ter regular intervals during the Monte Carlo simulation, 
all Nji{N}i — l)/2 pairs of the Nft rods are considered. If 
•s < Smax, where Smax IS a cut-off distance small compared 



// \ s ^ Smax, wnere Smax IS a cui-on aisiance smaii comparea 

(_[ J_dXij exp(— /3(7^[T/v] -I- J/bond)) ■ (4) to (Rg), rnj is incremented by one. Here, the indices i 



Here, C/bond is introduced as a constraint imposed on 
the bonds between neighboring hard spheres forming the 
membrane. It is zero if all bond lengths r^ are in the 
range a < 7'i < v3a and infinite otherwise. In our model, 
the only effect of adding inclusions to the membrane is 
to prohibit certain triangulations that do not adhere to 
the conditions described in the preceding section. In (|j), 
this can be accounted for by extending the definition of 
C^iond so that it diverges also if TJv cannot accommodate 
the Nr inclusions. 

Triangulations can be transformed into each other by 
means of the bond-flip algorithm and thus no vertex is 
different from the others. Consequently, the average en- 
ergy of a triangulation with, say, two inclusions only de- 
pends on their mutual distance s and angles 6i, 02, see 
Figure n^. This is the basic assumption of the present 
work. Strictly speaking, the distance s is the length of 
the shortest geodesic on the surface that connects the 
centers of the two inclusions. For s ^ ^, however, the 
Euclidean distance can be taken. The inclusion interac- 
tion is purely entropic and not a result of direct molecular 
interaction. However, an effective interaction potential 
^rr(s, ^1, ^2, k) can be defined as a potential of mean 
force (PMF) M by the relation 



$ 



RR 



$n 



-k^T In g 



(2) 



(5) 



where ^ma is the mutual avoidance potential of the rods 
and 5^^-' = g^^^[s,0i,92, k) is the pair distribution func- 
tion of the inclusions on the triangulation. 

The potential $rr was calculated analytically in g for 
the range (, ^ s ^ L. Surprisingly, the result depends 
on the sum of the angles di + 62 between the rods and 
the connecting vector. However, 9i + 62 is the same for 
the T-formation (Figure [l| c) and, for example, the case 
9i ~ O2 ^ 45° (Figure H b). For small distances s < L 
that are mainly considered in this article, the latter case 
is essentially a parallel side-by-side position, and thus 
parallel and perpendicular relative orientations (Figure 
|l| a,c) would be indistinguishable. For this reason, the 
angle difference 9 := \9i — ^2! is a better variable here 
for s K, L. On the other hand, one has to take into ac- 
count that for s > L, the parallel side-by-side position 
and e. g. the in-line positions (Figure |l| a,d) are then 
indistinguishable. 



and j are rounded to the closest integer i ~ 25 s/s„ 
and j ~ 25 9/{tt/2) (rods are head-tail symmetrical). For 
Smax, a length of 1.5L is used. The relationship between 



and g. 



(2) 



IS 



n. 



tj 



(A) 



(2) ^ 

^*^' ntot (A., 



(6) 



which is valid for s < Si„ax- Here, ritot is the total number 
of pairs of rods considered (some of which with s > Smax 
and thus ntot > J2i j ''^ij)^ ^-nd i^i) is the average surface 
area of a strip on the vesicle with s e [i, i -I- 1] x Smax/25, 
analogous to the area between two parallels on the globe. 




FIG. 4. If two idealized rods A and B have a center-center 
distance Pa Pb = s oi less than L, the usual accessible rota- 
tional angle of rod B around its center Pb of tt is reduced to 

■K - a ^ Pa+ Pb- 

The pair distribution function g"' is unity for a uni- 
form inclusion distribution on an exactly spherical vesi- 
cle. However, 

a) the vesicle shape deviates from the sphere especially 
for small k. 

b) Since inclusions are bound to vertices, "quantization" 
effects in the density can be expected if s becomes com- 
parable to (r) . 

c) In addition, a peak in g'^'^'^ for s k a (the smallest pos- 
sible distance) can be expected for the following reason: 
An inclusion straightens a row of vertices, thereby forc- 
ing its vicinity to remain closer to the perfect, flat and 
evenly spaced lattice. Trial moves that consist in plac- 
ing a rod parallel to an existing one might be less likely 
to violate a bond length condition because the lattice 
is more uniform there (closer to the regular hexagonal 
lattice). 

Effects b) and c) are typical lattice artifacts and there- 
fore unphysical. They are eliminated by the procedure 



outlined below. 

We are only interested in the variation in g^"^' caused 
by the increased probability of placing inclusions close 
to each other, for the membrane there is less corrugated 
normal to the plane than on average. In order to remove 
artifacts a-c, normalization runs are performed, accord- 
ing to the following prescription. The rod-leap MC step 
conforms in all details to the prescription described in 
the last section, except that before the bending ener- 
gies of the new and previous rod positions are compared, 
the new vertices are moved back to their original posi- 
tions. In other words, in the normalization run, each 
trial move is accepted if no bonds would be broken in 
moving the vertices to the positions Pn, but since ver- 
tices are never actually moved there, the lattice remains 
unchanged. This yields a normalization histogram riij 
and analogous to ([g) a pair distribution function g'^'. 
As a consequence, the moves in the normalization runs 
share all three types of artifacts a-c, but are insensitive 
to differences in bending energy. This allows to compute 
*rr(s,^) by 



$ 



RR, 



fcRT 



In 



(2) 
~(2) 

91 j 



In: 



nij "tot 



(7) 



Inclusions cannot overlap with others which implies a mu- 
tual avoidance condition of svnO > 2s/ L as shown in Fig- 
ure 0. Also, the center-center distance cannot be smaller 
than a hard sphere diameter, s > a. The mutual avoid- 
ance potential ^ma from (h) is infinite in these regions 
and zero elsewhere. 

Depletion effects ^,|l^ play no role in the present sim- 
ulation, since the lipids (free vertices) as the "small par- 
ticles" in the membrane network always retain their uni- 
form density. 



V. RESULTS 

A set of simulation runs is performed with different 
values for the bending stiffness coefficient k. For each 
simulation, a normalization run is performed accord- 
ing to the prescription outlined above. The simulation 
length is 10® MC-steps. Figure || shows — <i>R,R ^^(s, 0) 
for n/k-oT — 7.5. The rows in the back of the histogram 
{0 —^ 0) refer to rod pairs in parallel position; for small s, 
this must be a side-by-side arrangement (Figure |l| a,b). 
For s/L > 1, rod pairs in the in-line formation (Figure 
H d) also contribute to the average. The front rows in 
the diagram {9 — > 90°) refer to rods in the T- formation, 
which is slightly repulsive for very small s. This means 







that a close, perpendicular position reduces membrane 
fluctuations most efficiently. 




s/L 



FIG. 5. Plot of -$rr(s,6I) for iV = 1012 and 

K/fceT ~ 7.5. The plateau on the left side with the shape 
of the letter 'h' is the mutual avoidance area described in 
Sec. 



Ill; ^RR^s, 6) is not defined there and "I>ma diverges. The 



peak in the back corresponds to the interactive interaction of 
parallel side-by-side rod pairs, the slight dip in front to rods 
in the "T- formation" . 

In order to find out about the impact of many-body 
effects, the fraction of entries n^ stemming from isolated 
pairs of inclusions to those entries corresponding to clus- 
ters of three or more inclusions is calculated. An isolated 
pair is defined as a pair of inclusions with s < Smax in the 
absence of other rods which arc closer than Smax to either 
of the two. About 91% of the entries in riij stem from 
isolated pairs. Many-body effects are thus negligible in 
the present simulation as far as the attractive well of the 
potential is concerned. These results therefore confirm 
the conjecture stated at the end of Sec. Q. 

In Figure ^, —^kr{s,9 ~ 0) is shown for K/k^T = 
2.0, 7.5, and 11.0. The decay is very rapid with <I>rr 
vanishing approximately at s/L « 0.6 within the pre- 
cision of the simulation. At s/L — 0.75 — 3 (r) and 
s/L = 1.0 = 4 (r), remnants of the lattice periodicity in- 
duce small peaks in the graphs. The figure shows that 
<I>R,R is attractive for small distances s < L, but short- 
range. Netz 01 found <I'rr(s) to decay logarithmically. 
A fit with a logarithmic function of the type given in 0] 
is shown in Figure Hb ^ . In addition to the short-range 
attraction there might also be a long-range component. 
In p] fee r/ 128 was estimated for the magnitude of the 
leading term of long-range attraction in the effective po- 
tential between rodlikc inclusions. This value is less than 
the statistical error of our data in the limit of larger s 
(see Figure 1^). 



^This refers to inclusions with a quadratic term in the per- 
turbation Hamiltonian ItII. 
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FIG. 6. a) Plot of -$rr(s,6I = 0); the graph for 

K/ksT — 7.5 (+) corresponds to the shce from Figured along 
s at 6 = 0. The distance between two tics on the abscissa 
equals one mean bond length (r) (L = 4(r)). The vertical 
dashed line at s/L — a/L — 0.183 marks the smallest pos- 
sible bond length, b) as a) but for n/ksT = 11.0 together 
with a fit of ai ln(l — exp(— a2s/I/)), (ai — 7.5, a2 = 6.6) 
following the theory of Netz [7] . 
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FIG. 7. Plot of max(-$RR(s, 9 = 0)) for diflterent bend- 
ing moduli K. Inclusions on stifi'er membranes feel a more 
attractive potential, the relationship is roughly linear with 
max(-<l>RR(s,0)) « k/10. 

In other words, the effective attractive interaction be- 
tween a pair of inclusions becomes stronger the stifFcr the 
membrane of the vesicle is. This may seem contradictory 
because the stiffer the membrane the more suppressed 
are fluctuations of the membrane. However, we notice 
from the work of Netz [Q that short-range interactions 
between a pair of inclusions may be expected to depend 
logarithmically on the ratio s/^ for s/^ <C 1. A fit of 
the functional form for $rr(s) proposed by Netz to our 
data shows that the latter are consistent with a logarith- 
mic decay of the effective potential (see Figure gb) . The 
correlation length ^, on the other hand, increases with k 
as ^ ex exp(47rK/(3fcBr)) ||l^. Thus, combining the loga- 
rithmic decay of <i>RR with the latter expressions leads to 
a linear increase of — $rr with k for fixed s and T which 
is consistent with the plot in Figure R. 



VI. DISCUSSION AND CONCLUSIONS 



Consequently, our data do not permit us to comment 
on long-range effective interactions between inclusions. 
It is, however, noteworthy, that due to the computa- 
tional method of a perturbative treatment in terms of 
inverse powers of s, the results of MM cannot rule out 
the existence of short-range attractions consistent with 
the present findings. The results in |R], on the other 
hand, are exact and remain valid in the strong-coupling 
limit. 

The larger a, the more attractive becomes the inter- 
action. Figure [7| shows that this relationship is almost 
linear, with roughly max(— $rr(s, 9 = 0)) ~ k/10. 



The present article discusses Monte Carlo simula- 
tions for closed, triangulated membranes with mobile, 
rod-like inclusions as an extension of a widely used 
model |l6Hl9(| . Inclusions locally straighten the network, 
thereby quenching lateral fluctuations. The effective in- 
teraction potential is attractive and short-range (see Fig- 
ure Ig). For s — > 0, it is limited by the mutual avoidance 
condition (Figures ^ and ||) . 

The contribution of this simulation lies in the focus on 
small inclusion separations. We find an attractive inter- 
action potential of the order of fceT between rods that 
consist of five rigidly connected hard spheres, for bend- 



ing coefficients of n/k-oT of ttie order of 1 — 10 (Figure 
^). Furtliermore, the magnitude of the well depth of the 
potential grows almost Hnearly with k, as shown in Fig- 
ure 0. All of the results presented in this paper are valid 
in the dilute inclusion limit. The form of the function 
<i>rtR and its magnitude can be expected to change with 
higher densities. For very high inclusion densities, segre- 
gated phases might exist, one of them largely depleted of 
inclusions and another, dense phase with a nematic order 
due to the inclusion anisotropy. 

In conclusion, these results have implications for the 
formation of inclusion clusters. Such aggregation is fre- 
quently observed experimentally p6| and in simulation 
p7[ . A number of different driving mechanisms have 
been proposed, aside from those induced by direct inter- 
action or depletion forces p 10|. This includes sponta- 
neous curvatur e [ll|], l ateral tension y,g] , and conical in- 
clusion shapes ||3|.^lj|] . Fluctuation- induced interactions 
as considered in |^^ focus on the large-separation limit 
s/L :s> 1 and result in interaction potentials <I>rr ^ k^T. 
Although this is of great theoretical interest as it draws 
an analogy to the quantum-mechanical Casimir effect 
IgS],^, interaction energies far smaller than fceT are of 
little practical relevance in thermodynamic systems. The 
present work, in contrast, allows to gain insight into the 
strongly perturbed, short length-scale region, over dis- 
tances comparable to both the inclusion length and the 
lipid head size. This is the crucial length scale in chemical 
and biological applications. 
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